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SUMMARY 


An implicit finite-difference method is presented for obtaining 
steady-state solutions to the time-dependent, conservative Euler equa- 
tions for flows containing shocks. The method uses the two-point dif- 
ferencing approach of Keller with dissipation added at supersonic points 
via the retarded density concept. Application of the method to the one- 
dimensional nozzle flow equations for various combinations of subsonic 
and supersonic boundary conditions shows the method to be very efficient. 
Residuals are typically reduced to machine zero in approximately 35 time 
steps for 50 mesh points. It is shown that the scheme offers certain 
advantages over the more widely-used three-point schemes, especially in 
regard to application of boundary conditions. 

INTRODUCTION 

Over the past decade much progress has been achieved in computing 
mixed subsonic-supersonic flows using the potential equation. Recent 
results by Jameson (ref. 1) using a multi-grid algorithm demonstrate 
the efficiency now achievable in obtaining accurate solutions to the 
potential equation. In order to avoid the irrotational assumption 
inherent in the use of the potential equation, attention is now 
being directed toward developing efficient methods for solving the 
conservative Euler equations. Since the conservative Euler equations 



correctly model inviscid rotational flow and contain the Rankine- 
Hugoniot jump conditions, solutions with significant vorticity and/or 
strong shocks can be computed more accurately than those obtained 
using the potential equation. 

Stability and convergence of numerical solutions to the Euler equa- 
tions appear to be sensitive to the boundary conditions imposed. This 
is pointed out, for example, by Moretti (ref. 2) in regard to the cal- 
culations of Cline (ref. 3) for nozzle flows. This sensitivity is also 
seen in the numerical results of Yee, Beam, and Warming (ref. 4) who 
examined the effect of different boundary approximations on stability 
for implicit schemes which have as their basis a three-point central 
difference approximation for the Euler flux terms. Their study was 
motivated by the extra numerical boundary conditions, in addition to 
the physical ones, that are required by the three-point formulation in 
order to close the system of difference equations. This can be illus- 
trated using the one-dimensional Euler equations which consist of the 
conservation of mass, momentum, and energy. Written in nonconservat ion 
form, the characteristic slopes in the (x,t) plane are u, u i c. If 
one considers the case where the inflow and outflow conditions are sub- 
sonic, there will be two incoming characteristics at the inflow boundary 
and one at the outflow boundary, This implies that only three physical 
boundary conditions are necessary to solve the three partial-differential 
equations — two at the entrance and one at the exit. However, the use of 
a three-point difference representation for the flux terms requires six 
numerical boundary conditions because the difference equations can be 
applied only at the interior mesh points. Thus, either extra numerical 
boundary conditions or some other treatments at the boundaries are 
required. This necessity for extra conditions can be satisfied in 
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several ways; (1) Specifying all conditions at the boundaries, 

(2) additional differencing of selected conservation laws at the 
boundaries, (3) using spatial extrapolations or time-spatial extrap- 
olations at the boundaries, or (4) a combination of the above. These 
auxiliary conditions may be applied in an explicit or implicit manner. 
The paper by Yee, Beam, and Warming (ref. 4) examined the stability 
effects using all these possibilities. Every class of flow examined 
in reference 4 required extra numerical boundary conditions in addition 
to the physical ones. 

In order to avoid this need for extra numerical boundary condi- 
tions, we examine in this report the idea proposed by Keller (ref. 5) 
of writing any system of differential equations as a first-order 
differential system and then representing all derivatives by a two- 
point approximation; that is, (F ) = (F. - F. .,)/Ax. This concept 

Xji X 1*“X' 

i~ a 

was applied by Keller and Cebeci (ref. 6) to solve the boundary-layer 
equations, which contain terms like n . In order to write terms of 

w «/ 

this form as a first derivative, it was necessary to introduce a new 

variable, say s s u^, so that u^^ becomes s^. The resulting two- 

2 

point finite-difference equations were accurate to order Ay . Later, 
Wornom (ref. 7) provided a simple and efficient extension of Keller's 

4 

idea which is of order Ay accuracy for the boundary-layer equations. 
This extension still avoids the need for extra numerical boundary con- 
ditions required by three-point methods. 

The two-point differencing idea appears to be a very attractive 
choice for numerically solving natural first-order systems, such as 
the Euler equations, for several reasons. First, since the system 
contains only first derivatives, there is no need to introduce new 
dependent variables as was necessary for the boundary-layer equations 
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which contained second derivatives. Second, the two-point difference 
equations, when applied to solve the mass, momentum, and energy equa- 
tions without added dissipation terms, require only three boundary 
conditions to close the system; thus, for cases where three physical 
boundary conditions exist, the difficulties introduced at boundaries 
with a three-point method do not occur. The case of subsonic inflow 
and subsonic outflow is an example of a flow where three physical 
boundary conditions exist. However, there are also flow situations 
where the number of boundary conditions required to close the two- 
point system is larger than the number of physical boundary conditions 
available. These situations will also be addressed in this report. 

SYMBOLS 


^i’^i’^i 


^ 11 ’ ^ 12 ’ ^ 21' ^22 


^ 11 ’^ 12 ’^ 21 ’^22 


^11 

c 

d^l> d?2^> e<l) e^2> 
1 ’ ®i ’ ®i 

E 


cross-sectional nozzle area 
coefficients in block tri-diagonal system 
defined in equations (l9b) 
coefficients appearing in equations (l5) 
and (17) 

coefficients appearing in equations (15) 
and (17) 

coefficient appearing in equation (l5) 

2 

nofidimenslonal speed of sound, c * yP/p 

coefficients in solution algorithm;, defined in 
equations (20) 

nondimeusiGnal total energy times cross- 
sectional area 
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AD f(2) 

I 

M 

P 

R 

S 

t 

T 

u 

X 

p 

~p 

y 

Ap,Au, 

Ax 

At 

li 

Z . 


steady-state residuals; see equations (15) 
and (17) 

coefficients in boundary-condition equa- 
tion (26c) 

total number of mesh points 
Mach number 

nondimensional pressure times cross-sectional 
area 

residual norm; see equation (29) 

constant appearing in equation (24); related to 

entropy 

time coordinate 

nondimensional temperature (or enthalpy) 
nondimensional longitudinal velocity component 
axial space coordinate 

nondimensional density times cross-sectional 
area except in figures where it is the physical 
density 

retarded density times cross-sectional area; 

see equations (12) 

ratio of specific heats 

change with respect to time coordinate 

mesh spacing 

time step 

switch function defined in equation (12c) 
matrix defined by equation (20c) 
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Subscripts 

c 


i 

r 

total 

t 

X 

Superscripts 

n 

T 

(“) 

-1 


reference quantity in switch-function 
definition 

index counter in x-coordinate mesh 
reference quantity, nozzle entrance 
denotes upstream reservoir condition 
partial differentiation with respect to t 
partial differentiation with respect to x 

denotes time level 
denotes matrix transpose 
denotes retarded quantity 
denotes matrix inverse 
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GOVERNING EQUATIONS 


The governing equations are the quasi-one-dimensional , time- 
dependent Euler equations which define flow through a nozzle. These 
are written in conservation form as 

Pt + (Pu)jj = 0 (mass) (1) 

(pu)^ + (P + pu^)^ " I ® (momentum) (2) 


(P + E)uJ^ = 0 


(energy) 


(3) 


with 


P = (Y - 1) [e - I pu^] (gas law) (4) 

where p, P, and E are the density, pressure, and total energy times 

the nozzle cross-sectional area A(x). (See fig. 1.) p and u have 

2 

been nondimensionalized by p^, and u^ and P and E by 

where r denotes the entrance values. y is the ratio of specific 

heats. 

In this report only the steady-state solution is of interest, when 
it exists. Thus, the pressure is eliminated by integrating the steady- 
state energy equation taking the total enthalpy to be constant. The 
system of equations to be solved thus reduces to 

Pt + ( pu)x = 0 (mass) (5) 

( pu)^ + (P + pu^)j^ " f ^x (momentum) (6) 
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and 


P 


Y 


P 


(r 

\ total 



(7) 


where '^•total nondimensional reservoir stagnation temperature 

(or enthalpy). If equations (5) and (6) are written in quasilinear 
form, these equations will have two characteristics with slopes given 
by 


^ = (Y + l)u ± V(Y + l)^u^ - 4y(u^ - c^) 

dt 2 y (8) 


The directionality of the two characteristics is determined by the sign 
of equation (8) which is the same as the sign of u ± c where c is 
the local isentropic sound speed. The directionality of the charac- 
teristics at the boundary points will be used in a later section to 
determine how many boundary conditions may be applied at each boundary. 

FINITE-DIFFERENCE EQUATIONS 

The partial-differential equations (5) and (6) are replaced with 
the following two-point central-difference approximation: 


1 

2 






S3 


0 


(9) 


1 

2 


(PU), 


+ (Pu), 


i-1 


+ j(P + pu 


>1 - 


(p + 


PU^) 


-J, 


'Ax 



( 10 ) 


with 
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( 11 ) 


p. 

1 




Y 


Vtotal 


1 2 
2 


where Ax = constant = x. = (l-l)Ax, i = 1,2,... I. p is 

a retarded density used to introduce dissipation at supersonic points. 
This idea was developed independently by Hafez, South, and Murman 
(ref. 8) and Holst and Ballhaus (ref. 9) for obtaining solutions to 
the potential equation. The retarded density is given by 


P = p - ypj^Ax 


(12a) 


or 


Pi = Pi - - Pi_i) (12b) 

where 

= max |o,l - (M^/M^_^)^j (12c) 

is a reference Mach number and is the mid-cell value (taken 

as the average of and 

The retarded density method was chosen for two reasons. First, 
the solution algorithm can be written so that the extra mesh point 
introduced by the dissipation at supersonic points does not alter the 
solution algorithm; and, second, using the retarded density, dissipa- 
tion is only added at points where (supersonic or near- 

supersonic points). 

The coefficient y^^ given by equation (12c) is not unique for 
the Euler equations. Other forms for y^ were tried before selecting 
equation (12c). Equation (9) is linearized by Newton's method which 
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is achieved by letting 


P» = pH 1 + Ap 

(13a) 

u^ = u^“^ + Au 

(13b) 


and substituting these into equation (9). Neglecting terms of second 
order (e.g. , Ap • Au) this linearized equation can be written solely 
in terms of Ap and Au by eliminating the Ap terms using 


Ap^ = Apj^ “ ^i ^ ^'^Pi “ 

For simplicity the value of the switch function is frozen at the 

previous time step. It is multiplied by a density difference which 
is 0(Ax) except across discontinuities; ignoring the dependence of 
p on p and u in the linearization does not hinder convergence of 
the present algorithm. The desired approximation of the continuity 
equation is given by 


bij^APi + 


(15) 


where 


•’ll = |“i Mt] - ‘‘i 

_n-l/ 

bi2 = Pi /Ax 


a 
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(16a) 
(16 b) 
(16c) 


a 
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(16d) 
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'11 


- 1 -;;!/ 


Ax - 


2At 


n-1 

'^i-l 


(16e) 


f(l) = 
i 


vH-l 

(pu)i 


- (pu) 


T-l]h 


(16f) 


n 1 ri 1 

as + u“ ^P^- Inserting 


p" = p"-^ -H Ap 


To write the momentum equation in a similar form, we approximate 
the time term (pu) 

and u^ = u^~^ + Au and equation (11) into the momentum equation 
yields a linearized momentum equation with the same form as equa- 
tion (15), that is 


( 2 ) 

^21^^i ^22^^i ^21^^i-l ^22^^i-l ^i 


(17) 


The coefficients in equation (17) are left as an exercise to the reader. 

Note that equations (15) and (17) are compact in that only two 
points are used in subsonic regions and three points in supersonic 
regions (continuity equation only). In the next section, an algorithm 
is presented for solving the system of implicit difference equations, 
equations (15) and (17), assembled over all cells. 

SOLUTION ALGORITHM 

Equations (15) and (17) can be written as a block tridiagonal 
system which can be inverted using a block tridiagonal solver. Nor- 
mally, a tridiagonal system involves three mesh points where i+1, i, 
and i-1 would be the three diagonal terms. However, for the two- 
point system, the diagonal terms refer to different groupings of the 
dependent variables and not to separate mesh points. Block tridiagonal 
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solvers are applied in two recursive steps, one a forward sweep and the 
other a back substitution. 

The back-substitution step, used here is taken as: 

Apf = *^1^^ e^^^Au. (18a) 

i = I, I - 1, ... 2 

Au^_^ = d^^^ + e^^^Au^ (18b) 

This step is applied in a recursive manner for decreasing values of 
i (right-to-left sweep) and used the coefficients ®i^^» 

and which are computed from the forward step (left-to-right 

sweep ) . 

The block tridlagonal system inverted with this solver is obtained 
by writing equations (15) and (17) as 


^i(^Pl_2'^Pi-l^ ®i(^Pi 


AUi_j^) 


C^Au 


. = ffti), 

1 I 1 


( 2 ) 


T 


(19a) 


where 



Cii a^^^ 


^11 ^12 


^12 


^i = 


; = 


’ ^i " 


( 19 b ) 


0 &21 


_^21 ®" 22 _ 


? 22 , 



and (Ap. Ap . ^ )*^, ( Ap . , Au. ^ )*^, and Au. are the three diagonal 
terms. Other groupings of dependent variables are possible and, 
therefore, the algorithm given by equations (18) is not unique. 

The coefficients dj^\ and ep^ are obtained? by using 

equations (18) to eliminate the lower diagonal terms in equations (19), 
the coefficients are 

Cl ) ('2 -1 

®i ' ®i J ~ ^i ^i i = 2, 3, I - 1,1 (20a) 


12 



d 


( 1 ) 


d<2) 

1 


= zl\ 


(20b) 


where 


F, = 


,( 1 ) 


'11 


/ ,(1) . (1) (2)V ,(1) f(2) _ ,(1) 

^^^1-2 ®l-2^i-l/ ®'ll^i-l’ ®'21^i-l 


and 


= 


'11 

^21 


a + c e<^^e^2) + ^ ^d)' 

12 ®ll®i-2®i-l ^ll®i-l 

^ ( 1 ) 

^22 ^21®i-l 


(20c) 


These coefficients are computed in a recursive manner with increasing i 
(left-to-right sweep), the forward step. In order to start this sweep, 
the coefficients d^^\ ®0^^’ ^1^^’ must be known. They are 

obtained from the boundary conditions as discussed in the next section. 

After the coefficients are computed in the forward sweep, the back- 
substitution sweep is initiated using the outflow boundary condition to 
determine Auj. This is done by setting the terms resulting from the 
solution algorithm for Apj equal to those derived from the exit bound- 
ary condition. (See next section.) That is. 


Apj = dj^^ + Bj^^AUj “ Sj + hjAUj 

(algorithm) (boundary condition) 

Solving for AUj. we obtain 

Aui = [dj^> - gi]/[hi - ej^>] 


( 21 ) 


( 22 ) 


13 



BOUNDARY CONDITIONS 


Difference equations (15) and (17) are to be applied on the grid 
shown in figure 2 where i = 1 is the inflow boundary and i = I is 
the outflow boundary. The number of boundary conditions required to 
close the system of difference equations can be determined from the 
following relation: 

Total number of difference equations + number 

of boundary conditions = (number of dependent 

variables) • (number of mesh points) (23a) 

The two difference equations (mass, momentum) can be applied in 
all cells; therefore, the total number of difference equations is 
2(1 - 1). Thus, equation (23a) becomes 

2(1 - 1) + number of boundary conditions =2*1 (23b) 

and it appears that two boundary conditions are required to close the 
system of difference equations. However, if supersonic inflow exists, 
three boundary conditions will be required because to retard the den- 
sity at i = 1, the density at the "dummy" point i = 0 is required 
[^1 ~ *^1 ~ needed to compute using 

equation (12c). 

Next we examine whether the number of physical boundary conditions 
available are sufficient to satisfy the numerical boundary condi-* • 
tions required by the difference equations. 

Subsonic Inflow 

The characteristics for this case are shown in figure 3, where 
M = u/c is the Mach nximber and dt is a differential time. For 


14 



this case, there exists both an incoming and an outgoing characteris- 
tic. The single incoming characteristic can be replaced by an 
arbitrary boundary condition. Here, the boundary condition is taken 
to be 


(Pl/^l) /(p^/,4^)^ = S 


total 


(24) 


where S, . t is a constant based on reservoir conditions. Recall 
total 

that P and p are the physical pressure and density times the cross- 
sectional area. 


Subsonic Outflow 

The characteristics for this case are similar to the subsonic 
inflow. The incoming characteristic from downstream is replaced by 


Pt = P 
I exit 

where i = I is the location of the outflow boundary. 


(25) 


Supersonic Inflow 

The characteristics for this case are shown in figure 4. The two 
incoming characteristics can be replaced by two arbitrary boundary 
conditions. The entropy condition given in equation (24) is used as 
one of these boundary conditions and Pq is taken as the other. 

Supersonic Outflow 

The characteristics for this case are similar to the inflow case 
but are outgoing characteristics; and, therefore, no boundary condi- 
tions should be specified at this boundary. 
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Table I summarizes the number of physical boundary conditions 
available for four classes of flow as well as the number of boundary 
conditions required by the present two-point method. It can be seen 
from Table I that the two-point scheme is ideally suited for class' 1 
and 3 flows since the number of physical boundary conditions available 
equals the number of boundary conditions required to close the two- 
point system and these can be applied at the appropriate boundaries. 

For flows where the outflow is supersonic (classes 2 and 4), the two- 
point method requires one more boundary condition to close the system 
than there are physical boundary conditions available. The additional 
boundary condition needed for closure when the outflow is supersonic is 
taken to be the exit pressure. The consequence of this over-specifica- 
tion of the downstream boundary will be discussed in the section 
devoted to numerical results. 

In order to apply the boundary conditions to the difference equa- 
tions, they are written in the following forms; 


Apq = dQ^^ + Oq^^AUq (supersonic inflow only) 

APf = + e^^^AUj^ 

Apj ’ gj + hjAUj 


(26a) 

(26b) 

(26c) 


where d and e are superscripted so as to be consistent with the 
solution algorithm presented in the previous section.; Equation (26a) 
is applied as a Dirichlet boundary condition, = constant; thus 


d 


( 1 ) 

0 



0 
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In order to write boundary conditions (24) and (25) in the form 
of equations (26b) and (26c), we linearize the boundary conditions via 
Newton's method. This is accomplished by linearizing the boundary 
conditions (24) and (25) after eliminating the pressure with equa- 
. tion (11) and retaining only terms of order Ap and Au. Then the 
coefficients in equations (26b) and (26c) are: 


d 


( 1 ) 

1 





^total 



(28a) 


e 


( 1 ) 

1 



P^Ui/Dj 


(28b) 


Df ^total 


(28c) 


g 


I 



(28d) 




PlUj 



total 


1 2 \ 

- 2 


(28e) 


Recall that Auj is given by equation (22). All quantities in equa- 
tions (28) are evaluated at the previous time level. The superscript 
n-1 has been omitted for simplicity. 


NUMERICAL RESULTS 

The present method was applied to compute flow through the nozzles 
shown in figure 5. Table II summarizes the five flow conditions in- 
vestigated. These are the same cases investigated by Yee, Beam, 
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and Warming (ref. 4) using a three-point method. The initial condi- 
tions were obtained by linearly interpolating the velocity between the 
exact entrance and exit values and determining the density such that 
p. = (pu) . Linearly interpolating both density and velocity 

was also tried. The rate of convergence was unaffected by the choice 
of initial conditions. 

Subsonic Flow — No Shock (Case I) 

The convergence rate and Mach number profiles for the symmetric 
case are shown in figure 6. The residual shown in figure 6(a) and 
later figures is the maximum residual defined by 


R 


max 



max 



max 


(29) 


Cases 1 and II are ideal for the two-point scheme since the number of 
physical boundary conditions equals the number of boundary conditions 
required to close the system of difference equations and they are 
applied at the appropriate boundaries (cf. Table I, classes 1 and 3). 

This same case — when computed with a three-point method (ref. 4)— 
was the most sensitive (of all five cases) to the type of extra boundary 
conditions required by the three-point method. This is due to the 
throat being the only sonic point* Two types of extra numerical bound- 
ary conditions imposed on the three-point scheme produced converged 
solutions with shocks when a CFL (Courant, Fredricks, and LeWy) num- 
ber > 1 was used. A CFL^^^ = 1,000 calculation was reported with the 
three-point method with a profile given at 500 time steps. 

O 

The results shown in figure 6 were obtained with a CFL = 10 . 

max 

Accurate profiles were reached in 10 time steps. The reference Mach 
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number, M , used to compute the retarded density with equations (12) 
c 

was 1.0. When a value M =0.9 was used, the convergence rate was 

o 

unaffected but the local solutions at which - 0.9 were smoothed 
due to the added dissipation. 

Subsonic Flow — With Shock (Case II) 

The convergence rate, Mach number, and density profiles for this 

case are given in figure 7. The density Ptotal value of total 

density upstream of the shock. The density, p, shown in the figures is 

the physical value. This represents an inconsistency since previously 

it was defined as the physical value times the cross-sectional area. 

The value of M for this case and all other cases was 0.9. A value 

c 

of < 1 was necessary for transonic flows in order to introduce 
sufficient artificial viscosity at the sonic point to prevent an expan- 
sion shock from occurring there. The Mach number and density profiles 
are shown after 10 time steps. According to reference 4, the maximum 
CFL number for this case using a three-point method was 20; whereas, 

Q 

the present calculations were computed with a maximum CFL of 10 . How- 
ever, no noticeable change in the convergence rate was observed for 
CFL„„_ > 1,000 for all test calculations. On this grid (65 points), 
the residuals were reduced to machine zero in approximately 45 time 
steps with the largest residual always occurring at the shock. As in 
the other cases, plotting accuracy was achieved in 10 time steps. 

Subsonic Inflow — Supersonic Outflow (Case III) 

The convergence rate and density profile for this case are shown 
in figure 8. This case is of particular interest for the present method 
since a downstream boundary condition (pressure) is prescribed where 
physically none should be prescribed. The results presented in figure 8 
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demonstrate the method to be accurate and have fast convergence for 
this case. 

The results shown in figure 9 test the effect of applying a down- 
stream boundary condition by lowering the downstream pressure below 
the supersonic expansion value. Physically, any pressure below this 
value should have no effect on the solution in the nozzle. As seen 
from figure 9, the effect of this boundary condition is restricted to 
a small "boundary-layer" region near the outflow boundary. This 
boundary- layer region is caused by the artificial dissipation intro- 
duced by the retarded density. 

Supersonic Inflow — Subsonic Outflow (Case V) 

The convergence rate, Mach number, and density profiles for this 
case are shown in figure 10. As noted previously, in order to use the 
retarded density for flows where the inflow is supersonic, the density 
and Mach nvimber at the previous station must be specified. The values 
of Pq and Mq used were obtained from the exact solution. Again, 
according to reference 4, the maximiim experimental CFL number for the 
three-point method was 10. The experimental maximum CFL number for the 
present calculations was approximately 250. As with the other cases, 
plotting accuracy was achieved in 10 time steps and the residuals 
reduced to machine zero in approximately 35 time steps. 

Since this was the only test case for which the upper limit for 

a 

the CFL number was less than 10 (see Table III), more numerical ex- 
periments were carried out. It was found that virtually unlimited 
At could be used by starting with a moderate value of At (corres- 
ponding to, say, CFL = 250) and increasing rapidly. A few such 
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experiments showed convergence to machine zero could be obtained in 
less than 30 steps. 

Overall Comparison With Three-Point Method 

Table III shows a comparison between the experimental maximum 
stable CFL numbers for all five test cases for both the three-point 
method (ref. 4) and the present method. For this comparison, 50 mesh 
points were used with the following sampling of CFL numbers: 0.5, 1, 

10, 20, 10^, 250, 10^, and 10®. The maximum CFL shown for the three- 
point method is the maximum value reported for all choices for the 
extra numerical boundary conditions investigated for the three-point 
scheme. As shown in Table III, the present method permitted signifi- 
cantly higher CFL values to be used for three of the five cases. 

CONCLUSIONS 

An implicit finite-difference method has been presented for ob- 
taining steady-state solutions to the time-dependent quasi-one- 
dimensional Euler equations written in conservation form. The method 
uses a two-point difference scheme at subsonic points. At supersonic 
points, dissipation is added by the retarded density concept and the 
difference equations involve three points. However, the overall method 
retains the nice feature of general two-point methods as regards bound- 
ary conditions; that is, for subsonic inflow and outflow, the number of 
physical boundary conditions equals the number of boundary conditions 
required to solve the system of difference equations. In the present 
method, the outflow pressure is specified for all nozzle flows computed; 
thus, the method is very general in this respect. It was shown that 
this over-specification of the boundary conditions at a supersonic 
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outflow boundary was restricted to a thin "boundary-layer" region at 
the exit. 

In general, it is concluded that: 

(1) For subsonic inflow and outflow, with or without shocks, the 
present two-point scheme is preferred to a three-point scheme since it 
requires no extra numerical boundary conditions and these flows can be 

3 

computed at CFL numbers on the order of 10 larger than the three-point 
method. 

(2) For supersonic inflow, the present method requires boundary 
conditions at the first two entrance mesh points. 

(3) For three of the five test cases, the present two-point method 
permitted calculations with much larger CFL numbers than the three- 
point method. 
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TABLE I.- PHYSICAL AND NUMERICAL BOUNDARY CONDITIONS FOR VARIOUS FLOWS 


Class 

Description 


Boundary 

Conditions 


Physical 

Req'd. by Method 

Inflow 

Outflow 

Inflow 

Outflow 

1 

Subsonic inflow, subsonic outflow 

1 

1 

1 

1 

2 

Subsonfc inflow, supersonic outflow 

1 

0 

1 

1 

3 

Supersonic inflow, subsonic outflow 

2 

1 

2 

1 

4 

Supersonic inflow, supersonic outflow 

2 

0 

2 

1 




TABLE II.- SUMMARY OF FLOW CONDITIONS 


Case 

Nozzle Type 

Description 

Area Ratio 

I 

Convergent-divergent 

Subsonic inflow, subsonic outflow (no shock) 

2:1.16^ 

II 

Convergent-divergent 

Subsonic inflow, subsonic outflow (with shock) 

2. 5:1. 5 

III 

Convergent-divergent 

Subsonic inflow, supersonic outflow 

2:2 

IV 

Divergent 

Supersonic inflow, supersonic outflow 


V 

Divergent 

Supersonic inflow, subsonic outflow 



^Area ratio = 2:2 for the ssnnmetric case 




TABLE III.- EXPERIMENTAL MAXIMUM STABLE CFL NUMBER 


Case 

Present Method 

Three-Point Method 
(Ref. 4) 

I 

106 

10^ 

II 

10« 

20 

III 

10® 

10® 

IV 

10® 

10® 

V 

250^ 

10 


^Unlimited CFL number can be used with "gradual 
start." (See text regarding Case V.) 






Flow 



0 12 3 1-3 1-2 I-l I 

Figure 2.- Computational mesh. 
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i4(x) — ~ ^)/®j 

^(x) = " 5 )/(Xi-5)] 


(a) Convergent-divergent nozzle (ref. 10) ('Iq “ entrance area, 
= exit area, = throat area). 



i4(x) = 1.398 + 0.347 tanh (0.8x - 4) 


(b) Divergent nozzle (ref. 11). 


Figure 5.- Nozzle geometry. 
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Log,o(R) 



(a) Convergence history. 


Figure 6.- Case I: Subsonic inflow, subsonic outflow, no shock 

(CFLjjjax = 10® » 65-point grid). 
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2 


Figure 7,- 



Time Step 


(a) Convergence history. 


Case II: Subsonic inflow, subsonic outflow, with shock 

= 10^, 65-point grid). 
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btal 


Numerical solution ( 5 time steps) 
Exact solution 


□ 



X 

(b) Density. 
Figure 9.- Concluded. 
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